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Abstract 

We take a fresh look at the limits on the number of neutrino flavors derived from big bang 
nucleosynthesis. In particular, recent measurements of the 4 He abundance enable one to 
estimate the primordial 4 He mass fraction at Y p = 0.232 ± .003 (stat) ± .005 (syst). For a 
baryon to photon ratio, ip consistent with the other light elements, this leads to a best fit 
for the number of neutrino flavors N v < 3 (the precise number depends on rf) indicating 
a very strong upper limit to N v . Here, we derive new upper limits on N v , paying special 
attention to the fact that the best estimate may lie in an unphysical region (N v < 3 if all 
three neutrino flavors are light or massless; the lower bound to N u may even be as low as 2, 
if the small window for a v T mass is exploited). Our resulting upper limits therefore depend 
on whether N v > 2 or 3 is assumed. We also explore the sensitivity of our results to the 
adopted value of r/ and the assumed systematic errors in Y p . 


Among the strongest cosmological constraints on particle physics models are those de¬ 
rived from the consistency of standard big bang nucleosynthesis (SBBN) and, of particular 
importance are the constraints determined from the consistency of the predicted and ob¬ 
served primordial 4 He mass fraction, Y p . Most notable among these constraints is the limit 
on the number of neutrino flavors [|l] which can be translated into a host of other limits 
on particle physics properties. In the SBBN [Q, the abundances are primarily sensitive to 
a single parameter, the baryon-to-photon ratio, rj. Consistency between the predictions of 
SBBN and the observational determinations of the light element abundances restricts rj to a 
narrow range between 2.8 x ICC 10 and 4 x ICC 10 [||. In this range the calculated 4 He mass 
fraction lies in the range Y p = 0.240 — 0.245 |J, which is high when compared with the 
observationally inferred best primordial value 0 , 

Y p = 0.232 ± 0.003 ± 0.005 (1) 

where the errors are 1 a statistical and systematic errors respectively. Indeed, consistency 
within these small errors allows for very little room for any enhancement in primordial 4 He. 
This is the reason that SBBN leads to very tight limits on the number of neutrino flavors 

IBBB- 

Indeed, based on our estimate for the primordial mass fraction of 4 He, we deduced a best 
fit for the number of neutrino flavors |4J], 

N v = 2.17 ±0.27 ±0.42 (2) 

In H] we used a 2cr sta t + c S ys upper limit for F max in order to test the consistency of the 
SBBN. That procedure yields N v < 3.13, a value similar to, but somewhat weaker than the 
95 % CL limit N v < 3.04 from a Monte Carlo analysis ||. However, if there are indeed three 
massless neutrinos, then our best fit (N v = 2.17) is in fact in the unphysical regime. Limits 
derived from a probability distribution centered in an unphysical region of parameter space 
are known to give overly restrictive bounds.0 In this letter the limits on N u are rederived by 
renormalizing the probability distribution to the physical portion of parameter space. We 

1 We thank Paul Langacker for bringing this issue to our attention. 
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consider both N v > 3 and N u > 2 as criteria for the physical region. The limits on N v are of 
course weakened, even in the latter case (N v > 2) where though the mean value is physical, 
much of the probability distribution still lies in an unphysical region (Ay < 2). 

The best approach to a BBN bound on N v is to fit simultaneously both rj and N u to the 
inferred primordial abundances of all the light elements [;6|]. Here, instead, we simply explore 
these new limits to N u as a function of rj , as well as of the assumed systematic uncertainty. 
However, treatment of the systematic uncertainty can not always be done rigorously. We 
therefore study the implications of two different assumptions. We will assume either that they 
are Gaussian, in which case they will add in quadrature with the statistical errors or, that 
they are described by a top-hat distribution which we convolve with the Gaussian statistical 
errors (this is one of the cases being considered in 0). We also discuss the consequences of 
having several sources of systematic errors on our estimates of N u and Y max . 

As a prelude to our statistical analysis, it is useful to discuss first the relationship between 
N u and the 4 He abundance |lj . The 4 He abundance is primarily determined by the neutron- 
to-proton ratio just prior to nucleosynthesis or, more accurately, prior to the freeze-out of 
the weak interaction rates at a temperature Tf ~ 1 MeV. The final 4 He abundance is in fact 
quite sensitive to the freeze-out temperature which is determined by the competition between 
the weak interaction rates and the expansion rate of the Universe. As the uncertainty in 
the neutron mean-life is now very small, r n = 887 ± 2 s [7|], the helium abundance depends 
primarily on 7] which determines the onset of nucleosynthesis and N u , which can be used 
to characterize the expansion rate of the Universe at T > 1 MeV, as measured by the (time 
dependent) Hubble parameter H, 

H = / &**G„W. + lNr) \ 1/2 T , p) 

where iV* = 5.5 provides the contribution from electrons and photons. It is appropriate 
to emphasize that N v measures the contribution to the energy density at the epoch of 
nucleosynthesis of neutrinos (massive or massless) and any additional particles beyond those 
in the standard model. In the standard model, N u = 3 (unless the r neutrino has a mass in 
excess of about 0.1 MeV) [§. Limits on additional massless neutrino flavors with standard 
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model coupling to the Z gauge boson are very tightly constrained by LEP: N u < 3.04 (2 
cr) J7[]. The nucleosynthesis bounds however are most sensitive to any relativistic particle 
species present at that time. For example, an additional massless scalar degree of freedom 
such as the majoron, would contribute 4/7 to N v . 

The neutron-proton “freeze-out” is determined by the competition between the expansion 
rate (|j) and the weak interaction rate (proportional to T 5 ). A larger value for N u increases 
the expansion rate leading to earlier freeze-out at a higher neutron to proton ratio and, 
hence, yields a higher value for Y p . To understand the behavior of the calculated value of Y p 
as 7/io, r n and N v are varied, it is useful to fit the SBBN results in the form 


Y p — A + B In rjio + C(r n — 887) + D^Ny — 3) 


(4) 


As r/io increases from 2 to 10, B decreases from 0.014 to 0.009. For all ij, C = 2 x ICC 4 , 
confirming that the small uncertainty in the neutron lifetime contributes negligibly to the 
uncertainty in Y p . As N u varies from 2 to 4, D decreases from 0.014 to 0.012. From eq. (f|), it 
is clear that to bound N u will require that Y p and t/ 10 and their uncertainties be constrained. 

In Yang et al. 0 and Walker et al. we utilized solar system abundances of D and 
3 He to place an upper bound to the primordial abundance of D + 3 He and we derived 
a lower bound to 7/ 10 (> 2.8). This bound is conservative in the sense that it ignored 3 He 
production in low mass stars. Though models of galactic chemical evolution using primordial 
abundances near this lower bound on r] may show sufficient deuterium destruction m , when 
3 He production is included these models predict a large excess of 3 He 0] unless r] is large 
(£4)0. Indeed, a recent Monte Carlo analysis ||13|| found the larger value of ?/ 10 ~ 6.6±1.4 
as a best fit even when 3 He production was ignored. Thus, the challenge to BBN is clear 
since, even for ?/i 0 > 2.8 and r n > 885, we have Y p > 0.240 for N v — 3 . As was stressed 
in [Q, this is only consistent with the observational data if cr(syst) = 0.005 is added to the 
2cr stat upper bound of 0.238. 

Recently, there have been reports of observations of deuterium in quasar absorption 
systems 0 with a D/H abundance (which may be interpreted to be the primordial one) 
corresponding to a value of r/i 0 ~ 1.5. Although such a low value for r/ would be consistent 
with the observed 4 He abundance, it would exacerbate the problem due to the overproduction 
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of 3 He |I5| , We note that there is also the possibility that this observation can be 
interpreted as a H detection in which the absorber is displaced in velocity by 80 km s' 1 with 
respect to the quasar |15|]. Such an interpretation is likely, given recent reports [[17]] that a 
D/H measurement along a different line of sight may indicate a much lower value for D/H, 
corresponding to a much higher value for r/ 10 (> 6). If this latter is correct, it would appear 
to pose a conflict between the BBN predictions and the observational determinations of Y p . 
At this point, however, the data on D/H is too limited (and lacks consistency) to permit any 
firm conclusion. In our subsequent discussion of limits on N v , we will, where appropriate, 
give results for several values of ??io. 

Before returning to the problem of placing limits on N v , we first discuss our treatment 
of errors used in establishing a 2 a (or ~ 95% CL) upper limit. The chief problem is the 
treatment of systematic errors. In obtaining y max = 0.243 [f|], we simply added <r sys = 0.005 
to the 2a statistical upper limit. There are better alternatives. For example, it can be 
assumed (as is done by the particle data group 0) that systematic errors are Gaussian 
distributed and add in quadrature to the statistical errors. In this case, the 2 a upper limit 
to Y p is 0.244. Or, it may be assumed that the systematic errors are described by a top-hat 
distribution (ie. constant probability between ±cr sys and zero otherwise). The convolution 
of the Gaussian and top-hat distributions is just the difference of error functions given by 
(unnormalized) 


erf 


x T a S y S f-i 

v V^2^stat / 


-erf 


x — a 


sys 




\[2o s tat 


(5) 


where fi = 0.232 is the central value for Y p . For <r stat = 0.003 and a sys = 0.005, the 
distribution given in eq. ([5]) yields a 95 % CL upper limit of 0.240 which is somewhat more 
restrictive than the previous two estimates. 

It is very difficult to estimate the size of the systematic errors. Indeed, in exploring 
possible sources of such errors it is important to consider whether they may be correlated or 
uncorrelated among themselves |3| r |(§. It has recently been suggested that previous analyses 
have overlooked potential sources of error and, therefore, that the true systematic uncertainty 
may be considerably larger than our estimate of 0.005 ||T8| , |T9| . Although we do agree that 
the size of the systematic error remains uncertain, we consider that the recent analyses may 
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be overly simplified. To illustrate the problem, consider the contributions to the uncertainty 
from corrections for ionization and collisional excitation |3], ||. It is often assumed that there 
could be some neutral helium in the zone where hydrogen is fully ionized and, therefore, 
that the ionization correction will only increase the inferred helium abundance [[ 0 |, [19| . 
However, this overlooks the fact that very metal-poor stars will tend to be hotter than 
their higher metallicity counterparts, resulting in a harder ionizing radiation field in the 
very metal-poor HII regions. In such a situation, the Hell zone may extend beyond the 
HII zone and an appropriate ionization correction would reduce the inferred He abundance. 
Even worse, this possible effect would then correlate with the metallicity of the HII region, 
perhaps reducing Y at low metallicity and raising it for more metal-rich HII regions. The 
same considerations apply for the possible correction to account for collisional excitation in 
helium. This correction (decreasing the inferred helium abundance) will be larger for hotter 
(metal-poor) regions and smaller for cooler (metal-rich) ones. 

Thus, since some of the sources of systematic errors may be correlated (or anticorrelated) 
and, even their sign uncertain, a naive linear combination of several error estimates [p~8] , |H|j 
must lead to an overestimate of the true uncertainty. To explore the effect of such possible 
errors on our estimate of Y max , we consider some illustrative examples. For example, let us 
first assume that the systematic errors are Gaussian distributed. Then we simply add, in 
quadrature, the systematic error(s) to the statistical one. Thus, as noted above, with our 
choice of 0.005 for the systematic error and 0.003 for the statistical one, we find for the 95 
% CL upper limit on Y , T 95 = 0.244. For two uncorrelated systematic errors, each of size 
0.005, we would increase this to Y 95 = 0.247. If, instead, the systematic error were 0.010, 
we would find Y 95 = 0.253 ; note the relatively large difference if the systematic errors are 
added in quadrature or, linearly. In contrast, if we we treat the systematic errors as uniform 
(“top-hat”), the probability distribution for Y is given by eq. and for the three cases 
just considered we find Y 95 = 0.240, for a Gaussian statistical error of 0.003 and a top-hat 
systematic error of 0.005; 0.242, for a Gaussian statistical error of 0.003 and two top-hat 
systematic errors of 0.005 each (in this case one must convolve a second top hat with the 
distribution in eq. (|5|) which is a distribution given by a linear combination of four error 
functions and four Gaussians); and, 0.244 for a Gaussian statistical error of 0.003 and a 
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single top-hat systematic error of 0.010. These “2 cr” upper limits are considerably smaller 
than those from the double Gaussian approach and certainly much smaller than the upper 
bound obtained by adding the systematic errors linearly. 

We turn now to our computation of the upper limit on N v . To begin, we take Y p and its 
uncertainties from eq. (]4|) and compare to the central predicted value of Y p for each value 
of rjio and r n = 887 sec. to find the best (i.e., central) value for N u . For example, as 7710 
increases from 1 to 10, the central value of N v decreases from 3.6 to 1.1; the best values of 
N v corresponding to 7710 = 2.8 and 4.0, are 2.3 and 1.9 respectively. If all three neutrinos 
were in fact relativistic at the time of nucleosynthesis, then these latter values (and, indeed, 
all values for 7710 > 1.5) are unphysical. For the central case of 7710 = 3, even the 2 a statistical 
upper limit is unphysical. When the systematic uncertainty is included, the upper bound 
goes above 3.0, but this is most certainly an overly restrictive upper bound based on the data 
at hand. Here we will use the Bayesian approach described in [j7[]. Consider first, the double 
Gaussian treatment of errors (ie., both statistical and systematic errors as Gaussians). In 
the Gaussian distribution 

— (x — N u )‘ 2 

f{x) = e 2 ^ (6) 


a 2 = < 7 g tat + a 2 ys and N v is the (possibly unphysical) value determined for a given value of rj 10 . 
Normally, this distribution would be normalized by integrating over x from —00 to 00 and 
setting the result equal to unity. The 95 % CL limits on x would correspond to the limits 
of integration of the normalized distribution integrated to the value 0.95. Instead, here we 
allow the distribution to take non-zero values only over the physical region (N u > 3). The 
distribution is now normalized by integrating from 3 to 00 and the 95 % CL upper limit is 
the upper limit of integration X 95 so that the renormalized distribution yields a value 0.95 
when integrated from 3 to £ 95 , 


Is 95 f(x)dx 


= 0.95 


(7) 


JT f(x)dx 

The 95% CL upper limit on N v as a function of 7] w is shown in Figure 1 by the curve labeled 
N v > 3. For 77 10 > 2.8, N v < 3.61. This weaker bound would permit a massless scalar. 

As it is still possible that v T is massive and not relativistic at the time of nucleosynthesis, 
we also show in Figure 1 the result of the same calculation, when the probability distribution 
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is allowed to go down to N v = 2 (i.e. the lower limit of integration in eq. (0) is 2 rather than 
3). In this case, when 7710 > 2.8, N u < 3.18, a value closer to but still larger than the value 
of 3.13 found in the absence of renormalizing the distribution (even though the mean value 
is now physical, a large portion of the Gaussian is still unphysical). Now, a massless scalar 
is clearly excluded. It is interesting to note that, due to the renormalization procedure, the 
limit is actually weaker when three light or massless neutrinos are assumed. This is due to 
the renormalization procedure. 

In Figure 2, we show a similar plot with the systematic errors described by a top-hat 
distribution. The convolved Gaussian and top-hat is given by eq. ([5]) where now /j, corresponds 
to the mean value of N v and the errors are the propagated errors in N v rather than in Y p . 
The 95 % CL upper limits are found by using the convolved distribution (|5|) for /(x) in eq. 
d?]). The limits are now somewhat tighter at low 77 i 0 . In this case for r]i 0 > 2.8, N u < 2.91 
and 3.31 for N u > 2 and 3 respectively, excluding a massless scalar. 

Our results clearly depend on the choice of systematic error. In Figures 3 and 4, we 
show the sensitivity of our results to the value of (J sys . We plot the 95 % CL upper limit 
to N u as a function of <7 sys for N„ > 2 (solid curves) and N u > 3 (dashed curves) for three 
choices of 7710 = 1.5, 2.8 and 4.0. The curves are arranged such that the limits become more 
stringent (the upper limit on N v is lowered) as r]i 0 is increased. Figures 3 and 4 correspond 
to the double Gaussian and Gaussian/top-hat approach described above. Once again the 
latter provides a more stringent upper bound. As expected, as cr sys is increased the limits on 
N v become significantly weaker. 

In conclusion, we have shown that because the best fit values for the number of neutrino¬ 
like particle species is, or is close to being unphysical, the “true” upper limits on N v are 
significantly weaker than previous estimates. The upper limit on N v is quite dependent 
on whether one assumes N v > 2 or 3. Nevertheless, big bang nucleosynthesis can and 
does supply us with stringent constraint on N v and, more generally, on the speed-up of the 
expansion rate of the Universe. For the perhaps more likely choice of N v >3, and for 7710 in 
the range 2.8 - 4, the upper limit on N v ranges from 3.5 - 3.6 if systematic errors are treated 
as Gaussian and from 3.2 - 3.3 if systematic errors are treated as top-hats. We stress that 
these limits are only as good as the assumed errors (in particular the systematic errors); 
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the limit on N v is strongly dependent on <j sys . In the case of Gaussian systematic errors, 
as cr sys approaches 0.01, the bound on N v exceeds 4.0 (if N v > 3). If the systematic errors 
are treated as top-hats, then the bound on N v is less sensitive and only exceeds 4.0 when 
cr S y S > .02. Finally we note that even if the observational errors could eventually be greatly 
diminished, there remains a residual theoretical error < 0.001 in Y. Because this translates 
into an error of about 0.08 in N v , there will always remain a residual uncertainty of order 
0.08 in any determination of N v . Though we have argued here for somewhat weaker bounds 
on N v , they nevertheless still provide us with a strong means for limiting particle physics 
beyond the standard model. 
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Figure Captions 


Figure 1: The 95 % CL upper limit on N v as a function of 77 10 assuming that both statistical 
and systematic errors are Gaussian distributed. The two curves correspond to the 
physical condition that N v > 2 or > 3. In the former case it must be assumed 
that v r has a mass near its experimental upper limit. 

Figure 2: As in Figure 1, assuming that the systematic errors are described by a top-hat 
distribution. 

Figure 3: The 95 % CL upper limit on N v as a function of the systematic uncertainty in 
Y p . The solid (dashed) curves correspond to the condition that N u > 2 (3). Each 
of these two cases is shown for three choices of 7710 : 1.5, 2.8, and 4.0. The smaller 
values of 77 yield weaker upper limits . 

Figure 4: As in Figure 3, assuming that the systematic errors are described by a top-hat 
distribution. 
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